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Abstract 

We present calculations of the free energy, and hence the melting properties, of a simple tight- 
binding model for transition metals in the region of d-band filling near the middle of a d-series, 
the parameters of the model being designed to mimic molybdenum. The melting properties are 
calculated for pressures ranging from ambient to several Mbar. The model is intended to be 
the simplest possible tight-binding representation of the two basic parts of the energy: first, the 
pairwise repulsion due to Fermi exclusion; and second, the d-band bonding energy described in 
terms of an electronic density of states that depends on structure. In addition to the number of 
d-electrons, the model contains four parameters, which are adjusted to fit the pressure dependent 
d-band width and the zero-temperature pressure-volume relation of Mo. We show that the resulting 
model reproduces well the phonon dispersion relations of Mo in the body-centred-cubic structure, 
as well as the radial distribution function of the high-temperature solid and liquid given by earlier 
first-principles simulations. Our free-energy calculations start from the free energy of the liquid 
and solid phases of the purely repulsive pair-potential model, without d-band bonding. The free 
energy of the full tight-binding model is obtained from this by thermodynamic integration. The 
resulting melting properties of the model are quite close to those given by earlier first-principles 
work on Mo. An interpretation of these melting properties is provided by showing how they are 
related to those of the purely repulsive model. 
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I. INTRODUCTION 



Many years ago, a combination of experiments, first-principles calculations and simple 
models led to the comprehensive understanding of the low-temperature energetics of tran- 
sition metals that we have today- Much more recently, advances in experimental and first- 
principles techniques have started to open the possibility of achieving the same thing for 
the high-temperature phase diagrams, including melting curves, of transition metals over a 
wide range of pressures. However, the data obtained so far are fragmentary and sometimes 
conflicting,- and we believe there is now a clear need to develop simple models analogous to 
those used to interpret low temperature data. These models are needed in order to eluci- 
date the fundamental mechanisms that determine high-temperature phase diagrams, while 
providing a framework within which to interpret and unify experimental and first-principles 
data. We describe here how a simple parameterised tight-binding model can be used to 
calculate the high-temperature free energies of liquid and solid transition metals, and hence 
their melting properties, and we show how the model can help to interpret the available 
data. In the present work, we confine ourselves to the case of an approximately half-filled 
<i-band, focusing particularly on the interpretation of data for molybdenum. 

Shock measurements gave the first experimental information about melting curves at 
Mbar pressures, and data are available for several transition metals, including Fe, Mo, Ta 
and W.-^^ 1 ^ The thermodynamic states accessible in traditional shock experiments lie 
on a trajectory called the principal Hugoniot, which provides only a single point on the 
melting curve. On the other hand, major advances in static compression techniques, based 
on the diamond anvil cell (DAC), in principle allow entire melting curves and other phase 
boundaries to be mapped at pressures and temperatures up to ~ 200 GPa and ~ 4000 K. 
Melting data from static techniques have been reported for Fe, Mo, Ta, W, V and Y.-^^^^ 
There appear to be enormous differences between the melting curves of some transition 
metals from dynamic and static techniques, with the latter giving much lower melting slopes. 
The resulting differences of T m at Mbar pressures can be several thousand K. 

Melting curves from first-principles modelling began to appear over 10 years ago,— ^ and 
there are now several well established approaches, including the calculation of solid and liquid 
free energies, the "reference coexistence" method, and the explicit first-principles simulation 
of coexisting solid and liquid.— ii&iLi£iI2i22i2L22, For Fe, all three approaches have been used, 
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and the agreement between them is excellent.— 1 ^ Since DFT calculations are parameter- 
free, and reproduce very accurately key quantities such as cold compression curves, phonon 
frequencies, Hugoniot curves, and the zero-pressure melting temperatures of transition met- 
als, there is every reason to expect that their predictions of melting properties will also 
be reliable, and there is considerable evidence that this is the case. For transition metals 
for which static and dynamic measurements disagree seriously, first-principles calculations 
support the correctness of the dynamic measurements.— 1 ^ 

Molybdenum is one of the transition metals that have been intensively studied by DFT 
simulation, and it illustrates the recent controversies. Two independent sets of first-principles 
calculations^i^'S&SL^ agree rather closely with each other and support the high melting 
curve deduced from shock measurements, this curve rising far more steeply with pressure 
than the flat melting curve obtained from DAC data.— However, the shock measurements^ 
also indicate a solid-solid phase boundary, which may be the transition interpreted as melting 
in the DAC work.— A similar conflict between high shock and first-principles melting curves 
and a low DAC melting curve is also found in Ta,— an d it has been proposed that the 
transition seen in DAC may also be a solid-solid transition. We believe that simple models 
may help to resolve these controversies, by allowing the melting properties of transition 
metals to be related to the fundamental mechanisms that determine their energetics. 

Models for the energetics of transition metals are generally built on the principle that the 
total energy can be approximated as the sum of the electronic band energy and a repulsive 
pairwise interaction. The many different models that have been proposed differ mainly in 
their representation of the band energy. To explain the broad features of transition-metal 
energetics on a scale of several eV, including the roughly parabolic variations of cohesive 
energy, lattice parameter and bulk modulus with band filling, it suffices to assume a struc- 
tureless c?-band density of states (DOS), whose band width depends only on atomic volume 
(and chemical element).-^- The simplest total-energy model based on this idea consists of a 
sum of repulsive pair potentials plus a position-independent bonding term depending on the 
average atomic volume. We will refer to this as the REP+VOL model. More sophisticated 
types of total-energy models, including the closely related second- moment,— embedded- 
atom and Finnis-Sinclair models,— 1 ^ allow the second moment of the local DOS on each 
atom to depend on the distances to near neighbours. However, such models do not con- 
tain the physics needed to account for the well-known low-temperature structural sequence 
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that occurs through all the transition- metal series, from hexagonal close-packed (hep), to 
body-centred cubic (bec), to hep, and finally to face-centred cubic (fee). The energy dif- 
ferences of typically a few tenths of an eV between these structures are clearly essential 
for any discussion of phase diagrams, but they arise from the structure dependent form of 
the DOS. There are models that account for this by working with higher moments of the 
DOS than the second,- but a more straightforward approach is to express the total energy 
function directly in terms of a tight-binding (TB) model.— ^ In the present work, we use 
the simplest possible TB total-energy model, consisting of repulsive pair interactions plus 
the sum of single-electron energies calculated from a canonical <i-band TB model, without 
sp bands. We refer to this as the REP+TB model. With this simple model, we sacrifice the 
ability to describe the effect on the DOS of sp — d hybridisation, and the pressure dependent 
transfer of electrons between sp and d bands. We make this sacrifice in order to simplify 
the analysis. 

The principal question addressed in this paper is: What are the main parameters that 
determine the melting curves and other melting properties of transition metals, and what 
are the roles of these parameters? As part of this overall question, we would like to know 
at what level of detail we need to describe the (f-band bonding. In particular, do we need 
a detailed description of the structure-dependent electronic DOS in order to understand 
melting, or is a simpler model, such as REP+VOL, sufficient? In trying to answer these 
questions, our strategy will be to relate the melting properties of the REP+TB models to 
those of the pure REP model. 

Ultimately, we want to use parameterised tight-binding models to achieve a systematic 
overall understanding of the melting properties of the entire family of transition metals. 
However, even the simple models used here require rather extensive calculations to treat 
melting for a single metal, and for that reason we confine ourselves here to a narrow range 
of d-band filling in the region of half filling. We shall present a simple scheme for fixing the 
parameters of our model by fitting to zero-temperature first-principles data, and we shall see 
that, for the case of Mo treated here, we reproduce high-temperature first-principles results 
reasonably well. 

The remainder of the paper is organised as follows. In Section[TTl we present our REP+TB 
model for the total energy function, and we describe the scheme we use to fix the model 
parameters using information from T = K DFT calculations. In Section IHIl we present 
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a variety of tests of the model against DFT, both at T = K and for high-temperature 
solid and liquid Mo. The procedures used to calculate the free energies of the pure REP 
and REP+TB systems are described in Section HVl where we also report our results for the 
melting curves and the volume and entropy of melting. This is followed in Section [V] by 
an analysis of the relationships between the melting properties of the REP and REP+TB 
systems. Discussion and conclusions are in Section IVII . 

II. THE TIGHT-BINDING TOTAL-ENERGY MODEL 

The total energy U to t of our tight-binding (TB) model for a system of N atoms having 
position Ti is: 

tAotOi, r 2 , • • • r N ) = i ^ VwEpfaj) + U TB {ri, r 2 , . . . r N ) , (1) 

i¥=3 

where Vrep(?") is a repulsive pair potential and ry = [r» — r j | . In conventional TB treatments, 
the energy C/rB( r i> r 2 5 • • - r Jv) represents the sum of single-electron energies e n of occupied 
states, but here we include the effect of thermal excitation of electrons, so that Utb is 
actually a free energy, defined as: 

[/ T B = 2^/„e n -TS, (2) 

n 

where /„ is the Fermi-Dirac occupation number of energy eigenstate n at temperature T, 
and S is the electronic entropy, given by: 

S = 2k B J2 Ifn In f n + (1 - fn) ln(l - f n )] • (3) 

The factors of 2 in Eqns (T5]) and account for spin. The TB Hamiltonian used to calculate 
the e n is described next, and the repulsive pair potential is described after that. 

A. The canonical d-band tight-binding Hamiltonian 

Since we include only (^-electrons in our model, the Hamiltonian matrix elements 
(ia\H\j/3) (i, j label atoms, a, j3 label atomic orbitals) characterize hopping transitions 
of electrons between the d-orbitals xy, yz, zx, x 2 — y 2 , 3z 2 — r 2 on each atom. We employ 
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an orthogonal TB model, in which (ia\j/3) = 5ij5 a p. The dependence of the matrix elements 
on interatomic distance is taken to be exponential, so that: 

(ia\H\j/3) = Ga^r^Abexpi-rij/Rb) . (4) 

The factor G a p depends on the unit vector = (rj — rj)/rij in the direction from to ij, 
and it is well known that it can be expressed in terms of three basic matrix elements dder, 
dd7r and ddS. Here, we assume the canonical ratios^ ddcr : dd7r : ddS = —6 : 4 : —1. For 
convenience, and without loss of generality, we assume the diagonal elements (ia\H\ia) to 
be zero. In order to simplify the numerical simulations, we cut off the matrix elements so 
that they vanish beyond a distance -R cu t- The exponential is replaced by a cubic polynomial 
in the interval Ri < r < R cnt , the polynomial coefficients being chosen to ensure continuity 
of (ia\H\jf3) and its first derivative at R\ and -R C ut- For the Mo model developed here, we 
chose Ri = 4.7 A and R cut = 4.9 A. 

The TB density of states (DOS) n d (E), defined as: 

n 

is normalized so that J n d (E) dE = 10. Since the trace of (ia\H\jf3) is zero, the first moment 
[1^ of the DOS, defined as: 

ftf = J En d (E)dE j J n d (E)dE (6) 

(2) 

is zero. To fix the values of A b and Rb, we require that the second moment /i d of the DOS 
of our model, defined as: 

/4 2) = J E 2 n d (E)dE j J n d (E)dE , (7) 

should agree with the volume dependent <i-band second moment given by DFT. 

To apply this procedure to bcc Mo, we performed DFT calculations using the full- 
potential linearized augmented plane-wave method (FP-LAPW)— £2.>MM. as implemented in 
the WIEN2k code.— We used the Wu-Cohen 4 ^ form of generalized gradient approximation 
(GGA), which is known to perform well for transition metals.— 1 ^ Local orbitals are added 
to the standard LAPW basis in order to describe valence and semicore states. The technical 
parameters in the calculations were set as in Ref. 45| . The total and projected ci-channel 
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densities of states were obtained by using the modified tetrahedron method of Blochl et 
al.— , and for the projection we used an atomic sphere radius of typically 1.32 A. We found 
that Af, = 18.5745 eV and Rb = 0.8950 A give very good agreement with the DFT results 
for fj,^ at P = and 350 GPa (see Table I) and these values are used throughout this work. 

The quantity fi d is closely related to the width of the <i-band Wd, which is the difference 
between the lowest and highest energy levels, E d and E\ respectively, in the d-band DOS. 
In DFT calculations, the bottom of the (i-band E\ can be determined by direct inspection 
of the DOS, whereas E l d may be difficult to identify because of hybridization of states with 
different angular momenta (see Fig. [1]). Here we identify E\ with an abrupt drop in the 
projected g?-DOS at high energies followed by a smooth continuum. For Mo, we find that at 
equilibrium E\ and E d are —5.5 and 4.6 eV, respectively, while at a pressure of P = 350 GPa 
they are E\ = —10.8 and E d = 8.6 eV (see Fig. [Q. As shown in Table I, these values compare 
well with the tight-binding results obtained with the A b and Rb values quoted above. 

In order to reproduce the energy difference between the Fermi level and bottom of the 
(i-band and the form of the <i-DOS near Ep, we treat the number of d electrons as an 
adjustable parameter.— This is important, since many properties of transition metals are 
understood in terms of the form of the electronic DOS near the Fermi level (e.g. the relative 
stability of different structures, electronic specific heat, etc). For Mo, we find that Na = 4.3, 
rather than = 5.0, reproduces quite well the DFT results over a range of pressures (see 
Fig. CD). We use this value of Na, unless stated otherwise. 



B. The repulsive pair potential 

The pair potential VrepM is also assumed to have an exponential form: 

V REP (r) = A r exp(-r/R r ) . (8) 

The parameters A r and R r are chosen so as to reproduce as closely as possible the measured 
P-V curve of bcc Mo at low temperatures. This is essentially the same as fitting to DFT, 
since with the Wu-Cohen functional the DFT and experimental P-V curves are almost 
identical. The values A r = 3164.3454 eV and R r = 0.3350 A give excellent agreement with 
experimental data of Ref. jg], and DFT calculations (Fig. |2J), and we use them throughout 
this work. The same spatial cut-off distance and smoothing as used for the Hamiltonian 
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matrix elements is applied to the repulsive pair potential. 

III. SIMULATION TECHNIQUES AND TESTS OF THE MODEL 

A. Molecular dynamics simulation 

All the calculations on our TB model were performed with the OXON code,— 
using diagonalization of the Hamiltonian for each set of ionic positions. In the molecular 
dynamics (m.d.) simulations, we used the Verlet algorithm to integrate Newton's equations 
of motion, with a typical time step of 1.25 fs. The total force acting on each atom is 
the exact derivative of the total energy U tot with respect to its atomic position. Our m.d. 
simulations were performed in the canonical NVT ensemble, using Andersen's thermostat to 
avoid errors due to lack of ergodicity.— In using this thermostat, the atomic velocities were 
randomized by drawing them from a Maxwellian distribution every 0.2 ps. A typical m.d. 
run consisted of 2 ps for equilibration, followed by 10 ps for the calculation of averages. The 
m.d. simulations were performed on a 6 x 6 x 6 supercell containing N = 128 atoms, and 
T-point sampling was used to integrate over the first Brillouin zone. Pressure was obtained 
directly in each run using the virial formula. 

B. Tests of the model 

We have performed a series of zero and finite-temperature tests of our model in order 
to assess its accuracy compared with first-principles results and experimental data. In our 
first test, we evaluated the relative stability of the different crystal structures at different 
volumes. To this end, we computed the energy differences AE of the hep and fee structures 
with respect to the bec structure at equilibrium (Vq = 15.55 A 3 /atom) and a pressure 
of P = 350 GPa (V = 9.50 A 3 /atom) . At equilibrium, we found fcc-bcc and hep-bec 
differences of 0.40 and 0.46 eV with DFT, compared with 0.26 and 0.42 with TB. We thus 
reproduce correctly the relative zero-pressure stability of the different crystal structures, 
though we predict the fee phase to be appreciably more stable than hep. At P = 350 GPa, 
we found fcc-bcc and hep-bec differences of 0.30 and 0.32 eV with DFT, compared with 0.42 
and 0.30 eV with TB. At this pressure, we thus predict that the energy of the hep phase 
is lower than that of the fee phase, but bec is correctly predicted to be the most stable 
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structure. We note that at both pressures the agreement between the values of Ai^cp-bcc 
obtained with DFT and TB is very good. 

We have also computed the phonon frequencies of bcc Mo at the experimental equilib- 
rium volume Vq = 15.55 A 3 /atom .— Our calculations are based on the small-displacement 
method,— 1 ^ and we have used a supercell containing 64 atoms and 16 x 16 x 16 k-point 
grid over the first Brillouin zone. In Fig. [3] , we show our results together with experimental 
data from Ref. 6| and ab initio calculations from Ref. 25j for comparison; the agreement 
between the TB curves and the others is unexpectedly good, given the simplicity of our 
model. We note that the experimental phonon anomaly near the H point (1,0,0) is not well 
reproduced by either TB or DFT.— We have calculated the phonon frequencies also for 
the fee and hep structures of our TB model. We find that for = 4.3 and = 5.0 there 
are always imaginary frequencies, so that these structures are unstable, at least at T = 0. 
It worth noting that for slightly larger iV d values the fee and hep phases become stable at 
high pressures; for instance, for Na = 5.2 fee becomes stable at P ~ 400 GPa. 

The finite-temperature tests of our model include an analysis of the structure of the 
solid and liquid at different pressures. In Fig. H] , we plot the radial distribution function 
obtained from long (total simulation time ~ 10 ps) DFT and TB m.d. runs. The solid 
phase is simulated at T = 2000 K and P = 50 GPa , while the liquid is at T = 8250 K and 
P = 250 GPa. (These are states well below and above the melting curve of Mo given by 
first-principles calculations.}^ 1 ^ In both phases, the DFT and TB curves agree very well, 
the main difference being that TB gives interatomic distances slightly smaller than those 
from DFT simulations. 

In Fig. [5] , we show the electronic DOS of solid and liquid Mo obtained at the same 
thermodynamic conditions as for the radial distribution function. Although the DOS's 
obtained with TB and DFT are not identical, the corresponding band-widths and energy 
differences Ep — E b d are very similar, especially for the crystal. 

The main conclusion from all these tests is that, in spite of the formal simplicity of our 
TB model, it reproduces quite reliably many important properties of solid and liquid Mo. 
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IV. FREE ENERGY AND MELTING PROPERTIES OF THE MODEL 



Our overall strategy to obtain the melting properties of our model is based on the cal- 
culation of the Helmholtz free energy F tot (V,T) of the solid and liquid phases. To obtain 
Ftot{V, T) , we start from the Helmholtz free energy F REP (V, T) of the purely repulsive system 
described by the pair potential Vrep(^) , and use thermodynamic integration to determine 
the difference F tot (V,T) — F REP (V,T) at fixed (V, T). This thermodynamic integration is 
based on the general principle that for total-energy functions Uq and U\ , the difference of 
the corresponding free energies F (V,T) and Fi(V,T) at state point (V, T) is given by 

F x - F = [ (AU) X d\ , (9) 
Jo 

where AU = U\ — Uq and (•)> denotes the thermal average in the ensemble governed by the 
total energy function U\ = (1 — X)Uq + XUi . In practice, we use this type of thermody- 
namic integration to determine F tot (V,T) — F REP {V, T) at a set of (V,T) states, the value 
of F tot (V,T) at other states being obtained by integrating the relations P = —(dF tot /dV)T 
and E tot = (d(PF tot )/d(3)v, where (3 = l/k B T and E tot is the total internal energy of the 
system. The starting point of all the calculations is the free energy Frep(V,T) of the pure 
exponential system. Surprisingly, the thermodynamic properties of this system appear not 
to have been studied before, so we have performed our own calculations of -Frep(V, T), as 
described next. 



A. Free energy and phase diagram of the pure exponential model 

Details of the calculations of the free energy of the pure exponential model will be reported 
elsewhere, and here we give only a brief summary. Our values of -Frep(V, T) were obtained 
by thermodynamic integration (Eq. (j5J)), using as reference system the inverse-6 system 



interacting with pair potential VJ nv 6(V) = A/r 6 . We take the Helmho 



system from Ref. 



tz^free energy of this 
for the hep phase. 



571 ] for the liquid, bec and fee phases, and from Ref. 
The thermodynamic integration calculations were performed at a series of (V, T) points 
in which the free-energy difference -F REP — F inv6 was calculated by averaging Vrep — VJ nv6 
over long molecular dynamics runs in which U\ (see Eq. (Q) was varied continuously at a 
switching rate that guaranteed reversibility (that is, adiabatically). For the solid phase, we 
considered 12 volumes distributed uniformly over the interval 9.68 < V < 30.80 A 3 /atom, 
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and the temperature was set to T = 1000 K in all cases. For the liquid phase, 15 points 
within the same volume range as used for the solid and temperatures taken at intervals of 
1000 K from initial guessed melting temperatures up to 10000 K , were considered. We 
determined F tot (V, T) at the other state points by performing thermodynamic integration 
with respect to pressure and internal energy. 

Our calculated -Frep(V,T) values were cross-checked against simulations in which the 
liquid coexists with the bcc, fee or hep solid. These coexistence simulations were performed 



in the (JV, V, T) ensemble, and we used the techniques explained in Refs. [25] and [18j . 
Simulation boxes containing up to 10, 000 atoms were used in the calculations. For each 
pair of coexisting phases, the pressure dependence of the melting temperature T m (P) was 
fitted to the equation 



p y 
1 + T -1 



(10) 



b 

which resembles the so-called Simon equation^, but is adjusted to ensure that T m = at 
P = 0. Results for the bcc, hep and fee melting curves are shown in Fig. [6] . For the 
bcc melting curve, the values are accurately reproduced with parameters a = 564.6 K, 
b = 1.69 GPa and c = 0.5236. The volumes and enthalpies per atom of the coexisting 
solid and liquid were obtained from independent molecular dynamics simulations performed 
on supercells containing 1,000 atoms at (P,T) points on the melting curve. The melting 
volumes and entropy of fusion for the bcc melting curve are given in Table II. In fact, the 
melting curves obtained from the coexistence simulations were not perfectly consistent with 
the Helmholtz free energy results. We have searched carefully for the source of these errors, 
and we think it is possible that they may come from small imprecisions of the free energies 
of the inverse-6 system. To correct for these errors, we shifted the free energies of the 
bcc, hep and fee phases with respect to that of the liquid; the corrections depend solely on 
temperature, and are typically 10 — 20 meV/atom. We note that differences between energy 
shifts of all three crystal structures amount to less than 5 meV/atom, so therefore, total 
free energy differences between the bcc, fee, and hep phases, or equivalently their relative 
stability, are not affected appreciably by our corrections. 

As a further cross-check, we calculated the free energies of the bcc and fee phases by 
thermodynamic integration, starting from a different reference system. For this, we took 
a harmonically vibrating solid, with the harmonic force-constant matrix calculated for a 
system of particles interacting via the repulsive pure exponential pair potential at volume 
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V = 14.19 A 3 /atom (P ~ 204 GPa). These calculations are based on the small-displacement 
method.— 1 ^ For both the bcc and the fee phases, the small discrepancies between the free 
energies obtained from the inverse-6 and harmonic reference systems are typically 10 — 
20 meV/atom at temperatures near melting. 

B. Free energy and melting properties of the TB model 

Our thermodynamic integration calculations to determine the difference F to t — -Frep were 
performed by varying A adiabatically from to 1 over a time of 9 ps. Simulation boxes 
containing 128 atoms and T-point sampling over the first Brillouin zone were used. The 
quantity AU in these calculations (see Eq. ([9])) is the TB band free-energy [7 TB (see Eqs. (J2]) 
and ([3])). In order to reduce errors due to non-adiabaticity, we perform a complete cycle 
in which A goes from to 1 and back again, and to reduce statistical errors this whole 
cycle is repeated. The thermodynamic integral f dX (Utb)x is obtained as the average 
of the values in the four half-cycles. The typical standard deviation of these values is 
less than 10 meV/atom . An example of the (o 7 tb)a values obtained over a whole run at 

V = 15.55 A 3 /atom and temperature T = 5048 K is shown in Fig. [7J We note that (£/tb)a 
varies typically by ~ 0.5 eV/atom as A varies between and 1 . 

These thermodynamic integration calculations were performed at ten (V, T) states in 
each of the liquid and solid phases. For the solid phase, a temperature of T = 1000 K 
was chosen for all the volumes; volumes were drawn uniformly from the interval 9.68 - 
16.32 A 3 /atom. For the liquid phase, temperatures of typically 3000 K above the melting 
curve of the repulsive potential were chosen and the same set of volumes as for the solid 
was used. The value of F tot (V, T) at the other thermodynamic states was obtained by 
thermodynamic integration with respect to pressure and internal energy. 

In Fig. [8] , we report the melting line of our TB model for <i-band fillings Nd = 4.3 
and 5.0, obtained from the Helmholtz free energy calculations described above. We have 
considered different Nd values in order to assess the effect of this on the melting properties. 
Since our harmonic calculations showed that only the bcc structure is vibrationally stable, 
we will report only results for melting from the bcc structure. In practice, once the free 
energies F tot (V, T) of the liquid and solid phases are known, we have determined the melting 
pressure P m and volumes of the liquid and solid phases at each temperature by the Maxwell 
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double-tangent construction. The Simon formula T m = a (1 + P m /b) c was then used to fit 
our results and interpolate at any desired pressure. For Nd = 4.3(5.0), the values of the 
Simon parameters are a = 2865.9(1678.6) K, b = 118.3(35.1) GPa and c = 0.8530(0.6376). 
In Table III, we report results for the fractional change of volume AV/V S and entropy of 
fusion AS at points on the melting curves. 

In Fig. [HJ we also plot the melting curves of the pure exponential system and of Mo 
obtained from DFT calculations.— Although accurate reproduction of real-world data is not 
the main objective of this work, we note that our model (case Nd = 4.3) gives very good 
agreement with the P = melting temperatures for Mo of T m = 2883 K from experiment^, 
and T m = 2894 K from DFT simulations. 25 With increasing P, significant discrepancies 
between the TB (Nd = 4.3) and DFT melting curves appear. However, good agreement is 
partly recovered at high-P and high-T for Nd = 5.0. 

V. ANALYSIS OF MELTING RELATIONSHIPS 

We pointed out in the Introduction that the gross features of transition metal energetics 
at T = K can be understood on the basis of a model in which the structure of the 
electronic density of states DOS is ignored. This suggests that the simplest possible model 
for understanding the melting behaviour of transition metals is to add to the free energy of 
the pure exponential model Frep(V, T) a bonding term Ed(V) that depends only on volume 
and does not depend on temperature or on the phase of the system. To test this idea, we 
have carried out numerical calculations in which we have set Ed(V) equal to the bonding 
energy contribution to the total energy U to t of the bcc solid at zero temperature (A^ = 4.3). 
As expected, Ed(V) varies between —10 and —20 eV/atom over the volume range of interest. 
The resulting melting curve is shown in Fig. [9] . This very simple model necessarily shifts 
the melting curve upwards, and our results show that the computed melting temperatures 
are seriously overestimated, typically by around 50% . This result shows that there must 
be a significant dependence of the bonding energy on structure for given volume in the 
region of the melting curve. To illustrate this, we show in Fig. [TU] (Top) the bonding free 
energy AF = F tot — -Frep as a function of volume at T = 6000 K for the liquid and bcc 
solid. Remarkably, the difference between AF for liquid and solid is rather constant and 
has a value of ~ 0.2 eV/atom, AF being lower in the liquid. This means that the structure 
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dependence of the bonding stabilizes the liquid phase over the solid and therefore lowers the 
melting curve. 

It is interesting to ask whether AF is significantly influenced by the response of the struc- 
ture to the presence of the tight-binding energy. To answer this, we show in Fig. [TU] (Bottom) 
the quantity AF — (£/tb)rep where ({/tb)rep is the thermal average of U^b evaluated in the 
ensemble of the Vrep potential. The results show that AF — (?7tb)rep is quite significant 
in both phases. Moreover, the difference of this quantity for the liquid and solid indicates 
that the structure of the liquid responds significantly more than the solid to the presence of 
the TB energy. This effect contributes significantly to the lowering of the melting curve. 

VI. DISCUSSION AND CONCLUSIONS 

The present work is intended as a step towards developing an overall understanding 
of the phase diagrams of entire transition-metal series over a wide range of pressures and 
temperatures. At T = K, generalized phase diagrams (GPD) as a function of pressure F 
and atomic number Z can be computed by DFT, and we recently reported a phase diagram 
of this kind for the 4d series.— The construction of a complete GPD as a function of F, 
T and Z using DFT is too difficult at present, but we believe that it should be feasible 
using TB models of the kind described here. With this in mind, it is encouraging that our 
REP+TB model for Mo, parameterized using only T = K data, reproduces quite well the 
melting curve and properties of the high-T solid and liquid known from DFT. We have used 
the same REP+TB model, parameterized using the same scheme, for most of the other 4d 
metals, and we hope to report F-T phase diagrams for them in due course. We note that 
corresponding-states arguments will allow the free energies and melting data for the pure 
REP model reported here to be used to obtain the free energies of all these transition metals 
by thermodynamic integration. 

Since the properties of Mo over a wide range of F and T seem to be quite well described 
by our REP+TB model, it is natural to ask how the melting properties of the model are 
related to those of the pure REP model, consisting only of exponential repulsion. The 
melting temperature of REP goes to zero as F — > 0, so it is clear that the TB energy is 
crucial in determining the T m of transition metals at ambient F. However, we might expect 
the repulsion to be increasingly dominant at high F. There is an interesting connection 
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here with the melting properties of the Lennard- Jones (LJ) model for rare gases. It was 
recognised long ago^ that as P — > oo, the attractive r -6 potential of LJ has diminishing 
influence on the properties of the coexisting solid and liquid, so that the volume and entropy 
of fusion tend to those the soft-sphere repulsion r -12 model. As we have shown, the melting 
curves of our REP+TB model for Mo with Nj, = 4.3 and 5.0 do become close to that of pure 
REP at high P. Furthermore, the relative melting volumes AV/V S of REP+TB for both 
Nd values decrease steadily with increasing P, in a way that is consistent with convergence 
towards the melting volume of pure REP. However, this convergence is slow, since even at 
P ~ 400 GPa, AV/V S for REP is 1.1, while for REP+TB it is ~ 2.1 and 1.5 for N d = 4.3 
and 5.0 respectively. The entropy of fusion AS" is 0.74 k-Q for REP over the whole pressure 
range studied. For N d = 4.3, AS decreases steadily towards this value with increasing P, 
while for Nd = 5.0 it remains a little above this value for all P. 

In the Introduction, we asked what are the main parameters that determine the melting 
properties of transition metals. The success of our REP+TB model for Mo suggests that the 
two parameters A r and R r specifying the strength and range of the interatomic repulsion, 
the strength and range Af, and Rb of the TB matrix elements, and the number Nd of d- 
electrons, may be enough. (Firm conclusions must, of course, await TB calculations on 
other transition metals.) Our simulations show clearly that a description of the volume- 
dependent d-band width by itself is not enough. The very large d-bonding energy Ed(V) is 
described by the very simple REP+VOL model, but we have shown that this always raises 
the melting curve well above that of REP, and gives T m (P) predictions that agree poorly 
with the actual melting curves of REP+TB. The melting curves are substantially reduced 
below those of REP+VOL by the rather small shifts of relative free energies of solid and 
liquid included in the full REP+TB model. We have seen that a significant contribution to 
these shifts comes from the response of the system (particularly the liquid) to the presence 
of the TB energy. 

The present work may shed light on recent interpretations of the flat melting curves 
inferred from DAC measurements. It has been proposed that the directional bonding as- 
sociated with partially filled (i-bands may give rise to "preferred local structures" having 
icosahedral short-range order in the liquid phase.— ^ It was suggested that the formation 
of these local structures lowers the free energy of the liquid, and hence depresses T m . By 
contrast with simpler models, such as the embedded-atom model, the TB model we use 
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fully includes directional d-bonding, and our simulated liquid would presumably exhibit the 
effects of "preferred local structures", if they were present. The same can be said of the 
DFT simulations that have been reported on Mo. Nevertheless, both the present TB cal- 
culations and the earlier DFT simulations give much steeper melting curves than the DAC 
measurements, and this indicates that the full inclusion of directional d-bonding does not 
lead to low melting curves, in contradiction with the suggestions of Refs. 62.63] . 

Our TB model represents only the ci-band, and ignores the s — p band. This means that, 
although it mimics the pressure-dependent width of the (i-band and gives the main features 
of the DOS, it cannot reproduce the fine details, since it neglects hybridization of d-states 
with sp-states. It also means that the number of (^-electrons Nj, has to be treated as an 
adjustable parameter, and we do not include the dependence of this number on pressure or 
structure. It has been suggested 64 that the dependence of Na on structure might lead to 
the very flat melting curves inferred from DAC experiments. However, these ideas are not 
supported by DFT simulations, which fully include structure-dependent sp-d transfer, but 
nevertheless give melting curves that rise much more steeply than those from DAC. The fact 
that the present <i-band-only TB models give melting curves in reasonable agreement with 
DFT confirms that sp-d transfer is not expected to give flat melting curves. 

A possible resolution of the conflict between shock and first-principles melting curves on 
one side and DAC melting curves on the other side has emerged recently, at least for some 
transition metals.— DFT simulations of Mo have shown that, although bcc is the most 
stable structure at low T up to over 600 GPa, another structure, perhaps fee or hep, is likely 
to become more stable than bcc at much lower P and temperatures well below the melting 
curve. The suggestion is that the transition interpreted as melting in DAC experiments on 
Mo may actually be the transition between bcc and this other structure. 

Our main conclusions are as follows: A simple tight-binding model, parameterized using 
data for the volume-dependent <i-band width and the cold compression curve of Mo repro- 
duces reasonably well the melting curve and the properties of high-P/high-T solid and liquid 
Mo known from DFT simulations; the model allows us to analyse the physical mechanisms 
that determine the melting properties, and to assess suggested explanations for the anoma- 
lously low melting curves inferred from static compression experiments. We hope to report 
soon on TB calculations of melting properties across the whole 4d series. 
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w d 


(2) 


E F - E\ 


DFT 10.1(19.4) 


7.22(23.66) 


5.88(10.78) 


TB 10.8(19.4) 


7.29(21.25) 


5.85(10.86) 



TABLE I: Calculated ci-band width Wd = E l d — E b d , second moment /j, d and energy difference 
Ep — E b d from DFT and TB at P = GPa (P = 350 GPa in parentheses). Energies are in eV, and 
the number of d electrons is = 4.3. 
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P (GPa) 


T m (K) 


v t (A 3 ) f s (A 3 ) AF/y s (%) A5/fc B 


7.5 


800(100) 










58.7 


3100(100) 


20 74 


20 39 


1 73(5") 


73(2") 


88.6 


3985(100) 


18.45 


18.17 


1.54(5) 


0.74(2) 


141.1 


5200(100) 


16.09 


15.86 


1.49(5) 


0.74(2) 


204.5 


6400(100) 


1 4 38 


14 90 


1 98(5") 


73(9") 


269.0 


7450(100) 


13.23 


13.06 


1.32(5) 


0.75(2) 


333.5 


8450(100) 


12.36 


12.21 


1.20(5) 


0.74(2) 


409.5 


9450(100) 


11.57 


11.45 


1.07(5) 


0.74(2) 


491.5 


10450(100) 











TABLE II: Melting temperature function of pressure P, volumes per atom V[ and V s in 

coexisting liquid and solid, relative volume change AV/V S , and entropy of fusion AS of the pure 
exponential system for coexisting bcc solid and liquid. Estimated errors are given in parentheses. 
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P m (GPa) 

III \ "-*/ 


T (K) 


v, (A 3 ) 


v. (A 3 ) 

r s \ / 


AV/V. (%) 


AS/k B 




2000 


(16.99) 


(16.27) 


(4 37) 


(1 85) 


5 63(50 94) 


3000 


16 64(14 53) 


15 86(14 23) 


4 92(2 93) 


2 75(0 96) 


56.37(111.44) 


4000 


14.69(12.82) 


14.10(12.55) 


4.18(2.15) 


2.19(1.07) 


111.51(163.65) 


5000 


13.04(11.91) 


12.65(11.71) 


3.13(1.65) 


1.53(0.86) 


166.36(215.93) 


6000 


12.11(11.27) 


11.75(11.10) 


3.12(1.56) 


1.46(0.84) 


211.36(281.88) 


7000 


11.51(10.64) 


11.24(10.49) 


2.39(1.46) 


1.09(0.81) 


276.42(373.96) 


8000 


10.90(9.95) 


10.61(9.81) 


2.76(1.45) 


1.22(0.82) 


335.63(459.45) 


9000 


10.38(9.47) 


10.16(9.33) 


2.13(1.53) 


0.92(0.88) 



TABLE III: Melting function of temperature T, volumes per atom Vi and V, 

of coexisting liquid and solid, relative volume of fusion AV/V S , and entropy of fusion AS, for TB 
model at ci-band fillings Nd = 4.3 (5.0). 
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Figure Caption List 

FIG. [U : (i-component of the electronic density of states of bcc Mo calculated at T = K 
using DFT and TB at pressures of GPa (top panel) and 350 GPa (bottom panel). The 
Fermi energy is set to zero (vertical lines). 

FIG. [2] : Equation of state of bcc Mo obtained from the present TB model (solid line) 
and DFT (dashed line); experimental data (dots) from Ref. |6| are shown for comparison 

FIG. |3]: Phonon dispersion relations of bcc Mo calculated with the present tight-binding 
model (solid lines) and DFT (dashed lines) at the experimental equilibrium volume Vq = 
15.55 A -3 . Experimental data (dots) from Ref. 55] are shown for comparison. 

FIG. 1 : Radial distribution function of solid bcc Mo at T = 2000 K and P = 50 GPa 
from long DFT and TB m.d. runs. Bottom: Radial distribution function of liquid Mo at 
T = 8250 K and P = 250 GPa obtained from long DFT and TB m.d. runs. 

FIG. 0: (i-band electronic density of states calculated by DFT and TB m.d. simulation 
for bcc Mo at P = 50 GPa, T = 2000 K (top panel) and for liquid Mo at P = 250 GPa, 
T = 8250 K (bottom panel). 

FIG. [6] Phase diagram of the pure exponential model Vrep obtained from coexisting solid 
and liquid phase ( N, V, E ) simulations. The solid line in the figure corresponds to the 
bcc-liquid phase boundary while the dashed and dotted lines are the fcc-liquid and hcp- 
liquid ones, respectively. Dots simbolize points obtained directly from the phase coexistence 
simulations. 

FIG. [7] : Thermal average (£7tb)a of the tight-binding energy Utb as function of A in an 
adiabatic thermodynamic-integration calculation of the free energy difference F tot — -Frep 
between the REP+TB and REP systems. The plot shows {Utb)\ from a simulation in which 
A executes a double cycle O^l^O^l^O, the rate of variation \d\/dt\ being 1/9 ps _1 . 

FIG. E] : Melting curve of TB model at ci-band fillings N d = 4.3 (dashed line) and 5.0 
(dotted line). The melting curve of the pure exponential model and that of Mo from DFT 
simulations [25j] are show for comparison. 

FIG. M '■ Melting curve of the repulsive pure exponential potential Vrep (solid line), 
pure exponential potential plus a bonding energy term depending just on volume Vrep +E d 
(short-dashed line), and full tight-binding model U tot = Vrep + ^tb at N d = 4.3 (long-dashed 
line) and 5.0 (dotted line) . 
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FIG. [10] : Top: Free energy difference AF = F tot — -Frep of the total tight-binding 
model and repulsive pure exponential potential in the liquid and solid phases at temperature 
T = 6000 K and for N d = 4.3 . Bottom: Quantity AF — (?7tb)rep m the liquid and solid 
phases at temperature T = 6000 K and for iV d = 4.3 . 
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